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I.  INTRODUCTION 

The  need  for  the  analysis  presented  in  this  report  is  perhaps  best  described  by  Bill  Ormsby 
(Wll)  and  Jack  Carr  (ORCA  Computer  Inc.):  “An  acoustic,  passive  multistatic  network 
(A-PMSN)  for  situation  awareness  in  support  of  military  and  homeland  security  operations, 
such  as  harbor  defense,  is  needed.  The  A-PMSN  is  a  prototype  sound  source  location  system 
initially  focused  on  the  detection  of  small  boats  and  the  integration  of  those  detections 
into  existing  track  management  procedures  and  operational  picture  displays.  A-PMSN  is 
designed  to  operate  within  a  service-oriented  architecture  for  early  recognition  of  signihcant 
acoustic  detections  obtained  from  sensors  located  at  known  endpoints  of  a  baseline  using 
time  differences  of  arrivals  (TDOAs)  of  signals.  A  TDOA  is  dehned  to  be  the  difference  of 
times  of  arrival  of  a  signal  at  two  separated  receivers  (sensors). 

Using  sensors  on  more  than  one  independent  baseline  or  on  one  baseline  with  a  track 
constraint  allows  the  estimation  of  the  positions  of  the  sources  of  sounds  being  detected. 
These  positions  may  include  those  for  targets  of  interest  and  those  for  sources  of  sound  not 
of  interest.  The  TDOAs  will  be  produced  at  regular  time  intervals  and  converted  to  estimated 
positions.  These  positions  will  be  used  as  detections  and  processed  by  an  associated  cross 
correlation  algorithm  to  produce  tracks.  The  tracks  are  the  result  of  the  targets  moving 
from  hxed  point  to  hxed  point  over  the  span  of  the  data  collection.  In  particular,  targets 
whose  positions  fall  close  to  a  line  and/or  close  to  the  arc  of  a  circle  are  of  particular  interest 
and  form  the  basis  for  the  associated  algorithm  using  the  Hough  Transform  for  preliminary 
estimates  of  target  positions.” 

With  a  passive  multistatic  sensor  network  consisting  of  a  set  of  commercial  FM  or  acoustic 
transmitters,  receivers  (sensors),  and  a  data  processing  center  with  their  coordinate  locations, 
one  is  interested  in  locating  targets  moving  in  straight  lines  and/or  in  circular  or  possibly 
elliptical  arcs  within  a  given  surveillance  area.  An  appropriate  tool  to  aid  in  detecting 
preliminary  estimates  of  such  movements  is  the  Hough  Transform.  In  order  to  apply  the 
Hough  Transform  (HT),  N  input  points  (x,  y)  are  needed  that  are  generated  by  transmitter 
responses  from  a  target  to  receivers  (sensors)  using  TDOAs  [3] . 

Also,  two  methods  are  described  to  determine  the  (x,  y)  input  points.  First,  a  straight¬ 
forward  approach  dealing  directly  with  the  nonlinearity  of  the  problem  is  given,  followed  by 
a  particular  application  of  Ho-Chan’s  method  [4]. 

In  Section  H,  the  HT  is  dehned  and  discussed.  In  Section  HI  the  HT  is  used  to  identify 
targets  made  up  of  straight  lines  in  the  plane.  Sections  IV  and  V  describe  how  extensions  of 
the  HT  can  be  used  to  identify  targets  made  up  of  circular  and  elliptical  arcs,  respectively.  As 
noted  above,  the  required  input  is  a  set  of  possible  target  Cartesian  coordinates.  Numerical 
examples  are  given. 
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Section  VI  addresses  two  procedures  for  finding  the  possible  target  coordinates,  required 
by  the  HT  and  its  extensions,  using  TDOAs.  Three  numerical  examples  are  discussed. 

Fortran  95  software  has  been  developed  to  utilize  the  procedures  noted  above. 

II.  HOUGH  TRANSFORM 

Consider  a  straight  line  L  in  the  xy-plane  ( xyp  )  with  slope  m,  say 

y  =  m  (x  —  Xo  )  +  yo  ,  m  =  tan0,  (See  Figure  1),  (1) 

where  (xq,  yo)  refers  to  a  point  on  L.  Let  N  denote  the  straight  line  that  passes  through  the 

origin  and  is  perpendicular  to  L  at  a  unique  (xq,  yo)  as  shown  in  Figure  1.  Equation  (1),  can 
be  written  in  terms  of  the  parameters  that  describe  N,  namely  6  and  p.  The  result  is 

ysin0  +XCOS0  =  p,  — tt  <6<7r,  (2) 

where  6  denotes  the  angle  N  makes  with  the  positive  x-axis,  and  p  specifies  the  distance 

along  N  from  (0,  0)  to  (xo,  yo)  ■  If  in  (1),  (xo,  yo)  =  (0,  0),  then  p  =  0. 


Figure  1.  Normal  Form  of  the  Equation  for  a  Straight  Line 
Equation  (2)  is  referred  to  as  the  “normal  equation  for  L.”  It  is  obtained  from  (1)  by 


observing  that 

0  =  0±7r/2— )-m  =  tan0  =  —cos  9/ sin0  . 

(3) 

Hence  (1)  becomes 

y  =  -(cos6'/sin6')(x-Xo)  +  yo 

(4) 

y  sin  0  +  X  cos  6  =  cos  0  +  yo  sin  6  , 

(5) 

where 

P  = 

Xo  cos  0  +  yo  sin  0  =  9sz-\/ x^  +  yl  i  (See  Eigurel). 

(6) 
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It  may  happen  in  evaluating  p  from  (2)  that  p  <  0,  in  which  case  vr  is  added  to  6  if 
0  <  0,  or  TT  is  subtracted  from  0  if  0>  0,  thus  reversing  the  sign  of  p.  If  p  =  0  then  tt  is 
added  to  9  \{  9  <  0.  Hence,  it  is  assumed  that  p  is  always  nonnegative.  Note,  (2)  has  the 
advantage  over  (1)  since  no  problems  are  encountered  at  (j)  in  the  vicinity  of  ±7r/2. 

Consider  an  input  set  of  N  numerically  different  points  in  xyp  :  M  =  {Mn}  =  {(xn,  yn); 
n  =  1, . . .  ,N}.  Hereafter,  elements  of  M  will  be  designated  as  Mn  points.  Let  an  L  — line 
denote  a  straight  line  in  xyp  containing  at  least  three  Mn  points.  The  HT  transforms  L  to 
a  point  {9,  p)  in  the  0p-plane  {9pp),  where  the  point  is  taken  from  the  normal  form  of  the 
equation  for  L  as  described  above.  At  each  Mn,  a  spectrum  of  straight  lines  can  be  drawn 
by  varying  9  in  (2).  So,  each  such  Mn  generates  a  curve  in  0pp.  If  this  is  done  for  each  of  the 
N  points,  the  intersection  points  of  the  resulting  curves  in  0pp  designate  which  Mns  belong 
to  which  L  — lines,  thus  isolating  the  straight  lines  that  can  be  generated  from  the  N  given 
points.  Of  course,  one  would  look  for  three  or  more  intersecting  curves  at  a  point  in  0pp. 
Figures  2  and  3  illustrate  the  result  for  six  given  Mns.  Figure  3  exhibits  p<  0  and  applicable 
P>0. 


A 

O 


2 

O 


3 

O 


Figure  2.  Given  Mns,  (xn,  yn) ,  n  =  1, . . . ,  6,  in  xyp 


0— p  curves 
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Thus,  the  HT  can  be  used  to  isolate  sequences  of  points  that  he  on  straight  lines  embedded 
in  a  set  M  of  discrete  points  (xn,  yn),  n  =  1, . . . ,  N. 

The  graphical  approach,  of  using  the  curve  intersections  in  0pp  to  End  the  L  — lines,  i.e., 
lines  in  the  xyp  containing  more  than  two  Mns,  has  the  problem,  if  N  is  large,  that  the  curve 
intersections  in  0pp  are  obscure;  hence,  difficult  to  read.  In  such  cases,  if  one  wanted  the 
L— lines  containing  a  large  number  of  Mns,  it  could  be  difficult  to  separate  them  graphically. 
In  addition,  keep  in  mind  that  every  three  Mns  generate  an  intersection  of  two  curves  in  the 

0pp. 

Our  approach,  which  uses  no  graphing,  is  to  begin  by  finding  and  itemizing  the  Mns 
belonging  to  each  L  — line,  starting  with  the  hne(s)  with  the  largest  number  of  Mns  and 
listing  the  remaining  lines  in  decreasing  order  of  their  number  of  Mns.  More  precisely,  let 
Lk  denote  the  kth  L— line  obtained  from  the  HT  with  k  =  1,  2, . . .  ,  K,  and  let  Nk  denote  the 
total  number  of  Mns  contained  in  Lk.  Moreover,  the  K  Lks  are  ordered  in  decreasing  order 
of  Nk,  such  that  N1  >  N2  >,...,  NK  >  3. 

In  Figure  2,  the  Mns  are  indicated  by  small  numbered  circles.  Table  1  lists  the  coordinates 
of  these  points. 

Table  1.  Listing  of  Mns  of  Figure  2 

[xi,  X2,  X3,  X4,  X5,  Xg]  =  [  1,  5,  6,  3,  0,  7] 

[yi,  72,  ys,  74,  y5,  ye]  =  [1,7,1,4,7,10] 

Thus,  the  third  Mn  would  be  M3  =  (X3,  y3)  =  (6,  1). 

Using  Figures  2  and  3  as  an  example,  there  will  clearly  be  two  Lk  lines,  K  =  2.  So,  the 
output  would  be  as  shown  in  Table  2. 


Table  2.  Listing  Output  Two  L  — lines 


LI  ^ 

N1 

=  4, 

Mns  on  LI  — )■ 

M1,M2,M4,M6 

L2  -)■ 

N2 

=  3, 

Mns  on  L2  — )■ 

M3,  M4,  M5 

The  actual  mechanics  for  generating  software  to  determine  the  L  — lines,  Lk,  with  their 
Mns,  is  explained  below  using  the  example.  Every  two  Mns  determine  a  straight  line  with  an 
associated  unique  6  and  p  that  satisfy  (2);  then  a  6  and  p  are  determined  systematically  for 
all  possible  pairs  of  Mns.  Those  pairs  (at  least  two)  that  have  the  same  6  and  p  determine 
an  L  —  line  passing  through  those  Mns  (as  discussed  with  the  intersections  of  curves  in 
0pp  where  a  6  and  p  are  associated  with  each  intersection).  Beginning  the  process,  Mns 
one  and  two  would  generate  612  and  pi  2,  and  Mns  i  and  j  would  generate  j  and  p;  j 
with  i  <  j.  In  matrix  notation,  let  (0i,j )  and  (pij  )  denote  upper  triangular  matrices.  We 
use  a  four-quadrant  inverse  tangent  routine,  atan2(y,  x),  such  that  —tv  <  atan2(y,  x)  <  tv. 
Throughout,  the  subscripts  i,j  are  constrained  to 
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i<j,  i  =  l,2,...,(N-l),  j  =  2,3,...,N.  (7) 

With  e  discussed  at  the  end  of  this  section,  the  matrix  elements  are  generated  accordingly: 
If  |xi  -  Xj|  <  e  max(|xi|,  |xj|,  1)  and  |yi  -  yj|  >  e  max(|yi|,  |yj|,  l),then 

0,  if  Xi  >  0, 

TT,  if  Xi  <  0. 

Else,  if  |yi  — yjl  <emax(|yi|,  |yj|,  l)and|xi  — Xj|  >emax(|xi|,  |xj|,  l),then 

7r/2,  if  yi  >  0, 


Pij  =  |xi|  and  = 


Pi  i  =  lyil  and  d\  \  =  . 

'  -7r/2,  if  yi  <  0. 

Else,  if  (|xi  -  Xj|  >  e  max(|xi|,  |xj|,  1)  and  |yi  -  yj|  >  e  max(|yi|,  |yj|,  l)),then 

tan-i[(xi -Xj)/(yj -yi)],  i<j. 


^i,j  = 
A.j  = 


If  Ao  >  0, then 
else,  if  Pi  j  =0,  then 


0,  i>j, 

Xi  cos  0ij  +  yi  sin  0ij  ,  i<j, 

0,  i>j- 

PiJ  ~  Pi,]  1  ~  ^ij  ) 


(8) 


(9) 

(10) 

(11) 

(12) 


Aj  = 

=  A.j 

and  t 

r.  ^  j 

^  +  TT, 

if  Pi 
,  if 

A  IV 

(13) 

else,  if  Pi  j  <0,  then 

Aj  = 

-Aj 

and  1 

%]  =  < 

\  \] 

,  if  e 

,  if  e 

V  lA 

o 

(14) 

Referring  to  the  example  involving  Figures  2  and  3, 

we  obtain 

- 

-.5880 

1.5708  -.5880 

.1651 

-.5880\ 

0 

0 

.1651  -.5880 

1.5705 

;  -.5880 

0 

0 

0 

.7854 

.7854 

-.1107 

(15) 

0 

0 

0 

0 

.7854 

-.5880 

Vo 

0 

0 

0 

0 

1.9757/ 

/o  . 

2774 

1 

.2774 

1.1508 

.2774  \ 

0 

0 

6.0828 

.2774 

7 

.2774 

{Aj}  = 

0 

0 

0 

4.9497 

4.9497 

5.8527 

(16) 

0 

0 

0 

0 

4.9497 

.2774 

Vo 

0 

0 

0 

0 

6.4340/ 
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Note  that  =  tan^^[(0  —  7)/(10  —  7)]  =  —1.165904  with  Pg  g  =  ys  sin(— 1.165904)  = 
—6.43402.  Thus,  (14)  holds  and  to  insist  on  >  0,  ps^g  =  “Pgg  and  tv  is  added  to  ^g^g,  so 
that  ^5  g  =  05  g  +  TT  =  1.9757,  as  noted  in  (15).  The  software  which  generates  the  matrices 
in  (15)  and  (16)  is  the  subroutine  GUVM. 

The  software  designed  to  determine  the  Lk  lines  and  their  associated  Mns  is  Fortran  sub¬ 
routine  HSORTl  contained  in  hie  H:\CARR\MCARR.FOR  or  H:\CARR\MCARRl.FOR. 
The  reasoning  behind  the  software  is  explained  using  (8  —  16).  Starting  with  di^2,  one 
looks  along  the  hrst  row  for  0ij  =  0i^2  and  equality  in  the  corresponding  pij  =  pi,2- 
If  such  are  found  for  j  >  3,  then  one  need  look  no  further;  an  Lk  has  been  found.  From 
(12  —  13),  we  have  0i  2  =  0i,4  =  0i,g  =  —.5880  and  pi  2  =  Pi,4  =  Pi,g  =  .2774.  All 
the  subscripts  indicate  the  Mns  as  shown  in  Table  2.  There  is  no  need  to  examine  subse¬ 
quent  rows,  for  if  another  row  element  occurs  with  the  above  equalities,  its  Mn  will  already 
have  appeared  in  the  subscripts  obtained  from  the  original  row,  in  this  case  row  one.  For 
example  02,4  =  ^2,6  =  —-5880  and  p2,4  =  p2,g  =  .2774.  We  continue  with  searches  in  the 
hrst  row  of  (15),  but  no  more  matches  are  found.  In  the  third  row  we  End  two  matches: 
^3,4  =  ^3,5  =  -7854  and  p3^4  =  ps^g  =  4.9497  that  yield  the  Mns  given  in  second  row  of 
Table  2.  The  two  Lks  are  ordered  as  noted  above  and  listed  in  Table  2  as  LI,  L2.^ 

In  HSORTl,  numerical  equalities  for  the  elements  0i  j  and  pi  ^  are  dehned  in  terms  of 
an  input  epsilon,  e  >  0  (as  also  used  in  (8-11)),  as 


^ij  -  ®i.kl  £  max(l,  |«i,k|)£; 


(17) 


similarly  for  the  elements  pi j  and  pi^k  ,  numerical  equality  occurs  if 

|Pi,j  -  Pi,k|  <  max(l,  pij,  Pi,k)e.  (18) 

Two  references  with  an  extensive  discussion  of  HT  are  [6],  [7]. 


III.  TARGET  ID  USING  THE  HOUGH  TRANSFORM 

Two  more  complicated  applications  of  HT  are  discussed  here.  In  the  hrst  case,  we  consider 
16  Mn  points,  as  shown  in  Figure  4,  and  obtain  L-lines  using  HSORTl.  The  input  points, 
Mns,  are  also  given  in  Table  3. 


Table  3.  Listing  of  Mns  of  Figure  4 


[xi,  X2,  .  . 

•  ,Xlg] 

= 

[1 

5 

6 

3 

0 

7 

3 

3 

0 

5 

3 

7 

5 

0 

2 

5] 

[yi,  y2,  •  • 

■  Wie] 

= 

[1 

7 

1 

4 

7 

10 

1 

7 

0 

2 

3 

7 

5 

10 

0 

0] 

^ Speed  of  computation  was  necessary  for  HSORTl;  a  faster  intricate  program,  written  by  Russ  Gnoffo,  can  be  used  instead 
of  HSORTl  if  needed. 
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Figure  4.  Given  Mns,  (xn,  Yn) ,  n  =  1, . . .  ,  16,  in  xyp 


The  corresponding  curves  in  0pp  for  the  Mns  are  shown  in  Figured.  Note  the  relative 
difficulty  in  reading  the  intersection  points  in  0pp  . 


e 

Figure  5.  Transforms  of  the  Mns  in  Figured  to  0pp 

The  Lk  lines  determined  from  the  Mns  given  in  Table  3,  using  HSORTl,  are  given  in 
Table  4. 
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Table  4.  Listing  Output  Eleven  L  — lines  From  Table  3  Input 


LI  ^ 

N1  = 

5, 

Mns  on  LI  — ^ 

1,9,11,12,13 

L2  -)■ 

N2  = 

4, 

Mns  on  L2  — ^ 

1,2, 4,  6 

L3  ^ 

N3  = 

4, 

Mns  on  L3  — ^ 

4,7,8,11 

L4  — )■ 

N4  = 

4, 

Mns  on  L4  — )■ 

3,4,5,10 

L5  — )■ 

N5  = 

4, 

Mns  on  L5  — ^ 

2,10,13,16 

L6  ^ 

N6  = 

4, 

Mns  on  L6  — ^ 

2,5,8,12 

L7^ 

N7  = 

3, 

Mns  on  L7  — ^ 

1,3,7 

L8  ^ 

N8  = 

3, 

Mns  on  L8  — ^ 

4, 14, 16 

L9  ^ 

N9  = 

3, 

Mns  on  L9  — )■ 

8,13,14 

LIO  ^ 

NIO  = 

=  3, 

Mns  on  LIO  — )■ 

5,9,14 

Lll  ^ 

Nil  = 

=  3, 

Mns  on  Lll  — ^ 

9,15,16 

In  order  to  simplify  the  notation,  we  have  dropped  the  M  prefix  on  the  Mn  points,  so  that, 
for  example.  Ml,  M9,  Mil,  M12,  M13  for  LI  appear  as  noted  in  Table  4. 

The  second  application  of  this  section  involves  extracting  the  outline  of  a  house,  composed 
of  straight  lines,  emersed  in  a  maze  of  NR  =  350  random  points.  The  house  outline  is  made 
up  of  thirty  Mns  that  determine  five  L-lines.  The  points  are  given  by: 


<  y 

X 

.  y 


1.0, 

y  = 

=  1.0, 

1.2, 

1.4, 

1.6, 

1.8 

X  +  1, 

X  = 

=  1.0, 

1.2, 

1.4, 

1.6, 

1.8,  2.0 

5  —  X, 

X  = 

=  2.2, 

2.4, 

2.6, 

2.8 

3.0, 

y  = 

=  2.0, 

1.8, 

1.6, 

1.4, 

1.2 

1.0, 

X  = 

=  3.0, 

2.8, 

2.6, 

2.4, 

2.2,  2.0,  1.8,  1.6,  1.4,  1.2 

(19) 


The  NR  random  points  are  taken  from  a  Matlab  routine  based  on  normal  distributions. 
Therefore,  using  Matlab  code: 


mn  =  350; 
randn('seed',  1); 

<  X  =  .55  +  2.35  *  abs(randn(l,  mn));  (20) 

randn('seed',  2); 

y  =  .50  +  1.75  *  abs(randn(l,  mn)); 

Figure  6  shows  a  plot  of  the  30  Mns  and  those  random  points  in  the  interval  0.5  <  x  <  3.5. 
It  is  certainly  difficult  to  identify  the  house  outline  (30  Mns)  in  Figured.  Figure  7  shows  the 
same  points  with  the  Mn  points  darkened. 
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3.5 

3 

2.5 

2 

1 .5 

1 

0.5 

0.5  1  1.5  2  2.5  3  3.5 

Figure  6.  A  Plot  of  30  Mns  and  a  Subset  of  350  Random  Points 


3.5 

3 

2.5 

2 

1 .5 

1 

0.5 

0.5  1  1.5  2  2.5  3  3.5 

Figure?.  A  Plot  of  the  Identical  Data  of  Figure 6  with  the  30  Mns  Darkened 
The  input  array  (xj,  yi),  i  =  1, . . .  ,380  was  arranged  with  the  30  Mns  making  up  the 
hrst  30  points  of  the  array,  followed  by  the  NR  random  points,  but  the  results,  of  course, 
are  independent  of  the  order  in  which  the  input  points  are  given.  This  is  so  because  every 
possible  pair  of  points  in  the  input  array  is  treated. 


9 


NSWCDD/TR-08/136 


Calling  on  the  Fortran  subroutine  HSORTl  to  find  all  L-lines  containing  four  or  more 
Mns,  using  as  input  the  (x,  y)  of  (19)  and  (20),  the  output  is  given  in  Table  5. 


Table  5.  Listing  Output  from  Input  Given  by  (19)  and  (20) 


LI 

— )■ 

Nl= 

=11, 

Mns 

on 

LI 

— )■ 

1,  21,  22,  23,  24,  25,  26,  27,  28, 29,30 

L2 

— )■ 

N2= 

=6, 

Mns 

on 

L2 

— )■ 

11,12,13,14, 15,16 

L3 

— )■ 

N3= 

=6, 

Mns 

on 

L3 

— 

6,7,8,9,10,11 

L4 

— )■ 

N4= 

=6, 

Mns 

on 

L4 

— )■ 

1,2, 3, 4,  5,  6 

L5 

— )■ 

N5= 

=6, 

Mns 

on 

L5 

— )■ 

16,17, 18,19,20,21 

IV.  EXTENSION  OF  THE  HOUGH  TRANSFORM  TO  CIRCULAR  ARCS 

In  the  previous  section,  using  HT,  it  was  shown  that  Mns  to  determine  L-lines  could  be 
extracted  from  a  set  M  of  points  in  xyp.  The  principle  involved  is  simply  that  any  straight 
line  in  xyp  can  be  defined  by  two  parameters  {6,  p  in  the  HT  case),  so  that  examining  every 
pair  of  points  in  M  and  finding  those  pairs  with  the  same  9,  p  establishes  Mns  that  determine 
an  L-line. 

In  the  case  of  circular  arcs,  we  use  the  fact  that  a  circle  in  the  plane  is  uniquely  defined  by 
three  parameters.  Therefore,  to  determine  an  L  — arc,  i.e.,  an  arc  of  a  circle  (corresponding 
to  a  straight  line  in  HT  case),  one  would  look  at  all  possible  combinations  of  three  points 
from  M  and  search  for  those  particular  three-point  combinations  that  have  the  same  values 
for  each  of  the  three  parameters.  The  Fortran  software  that  accomplishes  this  is  contained 
in  file  DH.FOR.  Subprogram  DGXY  generates  M-sets  for  testing. 

We  use  for  the  three  parameters  defining  a  circle,  the  center  of  the  circle,  (a,  b),  and  its 


radius  r.  Thus,  the  equation  for  a  circle  in  xyp  is 

(x  -  a)2  +  (y  -  b)2  =  r^.  (21) 

The  objective  is  to  evaluate  a,  b,  and  r  using  three  distinct  points  from  M.  The  fact  that  a  and 
b  occur  nonlinearly  presents  no  problem.  Denote  the  three  points  by  (xl,  yl),  (x2,  y2),  (x3,  y3). 
Substitute  these  quantities  into  (21)  to  obtain 

(xl-a)2  +  (yl-b)2  =  r2,  (22) 

(x2-a)2  +  (y2-b)2  =  r2,  (23) 

(x3  -  a)2  +  (y3  -  b)2  =  r^.  (24) 

Then  subtracting  (22)  from  (23)  and  (22)  from  (24),  one  obtains,  after  cancelations: 

2  (x2  —  xl)  a  +  2  (y2  —  yl)  b  =  x2^  —  xl^  +  y2^  —  yl^,  (25) 

2  (x3  —  xl)  a  +  2  (y3  —  yl)  b  =  x3^  —  xl^  +  y3^  —  yl^.  (26) 
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Solving  the  last  two  equations  for  a  and  b  gives 

2a  =  [(x22+y22)(y3-yl)-(xl2+yl2)(y3-y2)-(x32  +  y32)(y2-yl)]/D,  (27) 
2b  =  [(x2^  +  y2^)(xl  —  x3)  +  (xl^  +  yl^)(x3  —  x2)  +  (x3^  +  y3^)(x2  —  xl)] /D,  (28) 

where  different  algebraic  expressions  for  D  are  given  by: 

(  (x2  -  xl) (y3  -  yl)  -  (x3  -  xl) (y2  -  yl)  ^  0, 

D=<  xl(y2  -  y3) +x2(y3  -  yl) +x3(yl  -  y2)  7^  0,  (29) 

y  yl(x3  —  x2)  +  y2(xl  —  x3)  +  y3(x2  —  xl)  7^  0. 

If  D  =  0,  then  either  the  three  input  points  lie  on  a  straight  line,  or  at  least  two  of  the 

points  are  the  same.  For  either  of  these  inputs,  no  circle  is  determined. 

The  remaining  unknown  r  is  obtained  from  any  one  of  (22)- (24),  using  a  and  b  obtained 
from  (27)  and  (28). 

The  software  that  reflects  this  analysis  is  the  subroutine  CIRl(xl,  yl,  x2,  y2,  x3,  y3,  a,  b,  r,  D). 
Hence,  to  extract  the  circular  arcs  from  M,  the  subroutine  CIRSOR  generates  three  arrays 
U(i,j,k),  V(i,j,k),  and  R(i,j,k)  that  are  functions  of  triplets  of  all  possible  combinations  of 
three  points  from  M,  without  repetition  (of  which  there  are  N!/[3!(N-3)!]).  So  that  for  a 
hxed  i,j,k  the  element  of  each  array  is  constructed  by  CIRSOR  as: 


U(i,  j,  k) 

=  a. 

i  =  1,. 

..,N-2, 

j  =  i  +  l,- 

.,N-1, 

k  =  j  +  1, . 

.,N, 

V(i,  j,  k) 

=  b. 

i  =  1, . 

..,N-2, 

j  =  i  +  l,- 

.,N-1, 

k  =  j  +  1, . 

.,N, 

R(i,  j,  k) 

= 

i  =  1, . 

..,N-2, 

j  =  i  +  1,- 

.,N-1, 

k  =  j  +  1, . 

.,N, 

where  CIRSOR  calls  CIRl  to  hud  the  elements  a,  b,  r^  of  U,  V,  R,  respectively,  corresponding 
to  each  triplet  of  input.  If  D  =  0,  r  is  set  to  zero. 

The  similarity  to  (10)  and  (11)  should  be  clear. 

As  a  special  case,  hie  PTCIR.FOR  contains  Fortran  software  that  Ends  L  — lines  each 
with  an  attached  L— arc  at  an  endpoint.  As  an  example,  we  consider  M  made  up  of  a  set  of 
NR=400  random  points  from  normal  distributions,  a  horizontal  line  L  with  NL=19  points, 
and  a  circular  arc,  L  — arc,  with  NC=10  points,  that  begins  at  the  right  end-point  of  L, 
(y  =  3,  1  <  X  <  3)  and  forms  an  arc  of  the  circle:  (x  —  1)^  +  (y  —  4)^  =  5.  Hence  M 
contains  a  total  of  N=429  points  in  xyp,  where  the  Fortran  95  code  that  follows  indicates 
the  generation  of  the  NR=400  random  points,  followed  by  the  NL=19  straight  line  points, 
and  ends  by  generating  the  ten  L— arc  points.  The  code  that  generates  these  points  is  carried 
out  by  the  subroutine  GXYC,  which  follows. 
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DIMENSION  X1(429),Y1(429) 

NR  =  400 
NL  =  19 
NC  =  10 

IR  =  1  !Seed  for  random  number  generator. 

IS  =  2 

CALLRN0R(IR,X1,NR,IER)  ISubroutine  RNOR  taken  from  [5], 

CALL  RN0R(IS,Y1,NR,IES)  IRandom  coordinates. 

DO  5  I  =  1,NR 

X1(I)  =  3.5*ABS(X1(I)) 

5  Y1(I)  =  4.  +  Y1(I) 

AO  =  1. 

BOl  =  3. 

DT  =  (BOl  -  AO)/NL 
NRl  =  NR  +  1 
Xl(NRl)  =  AO  +  DT 
Yl(NRl)  =  BOl 
DO  10  I  =  NRl,  NR+NL  -  1 
X1(I+1)  =  X1(I)  +  DT 
10  Y1(I+1)  =  Y1(I) 

PI  =  3.14159265359 
NRC  =  NR  +  NL  +  1 
NCI  =  (NC  -  1)*2 
PIN  =  PI/NCI 
BO  =  4. 

RO  =  ^ 

THETA  =  ATAN2((B01-B0),(B01-A0)) 

DO  15  I  =  NRC,N 

THETA  =  THETA  +  PIN 
X1(I)  =  AO  +  RO*COS(THETA) 

15  Y1(I)  =  BO  +  RO*SIN(THETA) 

Figure  8  shows  some  of  the  random  points  and  the  points  belonging  to  the  straight  line,  but 
the  points  belonging  to  the  arc,  although  included,  are  not  evident.  Figure  9  has  the  same 
points  as  Figure  8  with  the  straight  line  and  circular  arc  points  enhanced. 
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Figures.  A  Plot  of  Random  Data  and  Some  Points  Forming  a  Straight  Line 

with  an  Attached  Circular  Arc 
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Figure  9.  Same  Plot  as  in  Fignre  8  with  the  Straight  Line  and  Circular  Arc  Points  Enhanced 

The  object  of  this  simulation  is  now  to  extract  the  Mn  points  belonging  to  the  L-line 
using  subroutine  HSORTl  and  the  attached  L  — arc  using  subroutine  HSORT2.  The  results 
are  given  in  Table  6. 
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Table  6.  Listing  Output  Using  HSORTl  and  HSORT2 

***DOUBLE  PRECISION  PTCIRC.FOR  ***  JM=  0,  MZ=  4 

NR  =  400,  NL  =  19,  NC  =  10,  N  =  429,  AO=l,  BO=4,  RO=2.23607 

1  19  (The  N  points  of  M  are  input  for  HSORTl.) 

401  402  403  404  405  406  407  408  409  410  411  412  413  414  415  416  417  418  419 

IX  419 

1  10  (Only  the  N  points  of  M  and  attached  point  419  are  input  for  HSORT2.) 

420  421  422  423  424  425  426  427  428  429 

2  6 

123  154  189  239  291  298 

3  6 

53  180  190  233  306  329 

In  the  hrst  line  of  the  table,  JM  is  not  relevant;  MZ=4  indicates  that  only  output  with 
at  least  four  Mu’s  will  be  given.  The  second  line  specihes  the  input  for  generating  the  test 
data.  The  third  line  notes  the  hrst  line  of  output,  LI,  using  HSORTl,  and  indicates  19 
points.  The  fourth  line  lists  the  19  points.  The  sixth  line  indicates  the  point  at  which  the 
L-line  and  L— arc  are  attached.  The  seventh  line  gives  the  hrst  line  of  output,  with  10  points, 
using  HSORT2,  referring  to  the  L— arc  .  The  eighth  line  specihcally  gives  the  10  points.  The 
last  four  lines  refer  to  L  — arcs  contained  in  the  distribution  of  random  points  and  attached 
to  point  419. 

Note,  if  (xe,ye)  denotes  the  last  point  of  L,  that  all  triplets  of  M  chosen  include  (xe,ye). 
Hence,  only  every  pair  of  points  of  M,  without  repetition,  make  up  the  input  together  with 
(xe,ye),  rather  than  an  independent  triplet  of  (x,y)  points.  The  software  returns  results  more 
quickly  under  these  conditions. 

In  the  hnal  example  of  this  section,  it  is  required  to  hud  points  belonging  to  three  circular 
arcs,  three  L— arcs  ,  embedded  in  a  set  M  also  containing  random  points.  This  dihers  from  the 
previous  example  where  the  L  — arc  was  attached  to  the  L-line.  Thus,  here  the  three  points 
required  to  establish  an  L  — arc  are  independent  so  that  the  size  of  the  problem  increases 
signihcantly  from  an  N-squared  problem,  as  in  the  previous  example,  to  an  N-cubed  problem. 
We  need  to  consider  all  possible  combinations  of  the  points  of  M  taking  three  points  at  a 
time  without  repetitions.  The  total  number  of  such  combinations  is  NT=N(N-l)(N-2)/6 
rather  than  N(N-l)/2  as  in  the  previous  example  where  all  possible  pairs  were  considered. 

File  DH.FOR  contains  subroutine  DGXY  that  constructs  the  simulation  by  specifying 
the  three  L  — arcs ;  L  —  ai,  a  =  aoi,  b  =  boi,  r  =  roi,  i  =  1,  2, 3.  Thus 
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L 


L 


al  :  X  =  aol  +  rol  *  cos  6, 
y  =  bol  +  rol  *  sin  9, 

a2  :  X  =  ao2  +  ro2  *  cos  6, 
y  =  bo2  +  ro2  *  sin  6, 


aol  =  1,  bol  =  4,  rol  =  4 

30  <  0  <  100  (degrees),  11  plotted  points. 

ao2  =  3,  bo2  =  4,  ro2  =  3 

30  <  9  <  150  (degrees),  21  plotted  points. 


L  —  a3  :  x  =  ao3  +  ro3  *  cos  9, 
y  =  bo3  +  ro3  *  sin  9, 


ao3  =  —1,  bo3  =  4,  ro3  =  3.5 

-60  <  0  <  60  (degrees),  16  plotted  points. 


(31) 

(32) 

(33) 


Figures  10  and  11  show  the  points  of  M  within  the  x  and  y  intervals  of  the  graphs  where 
N  =  160  and  NT  =  N  (N  -  1)  (N  -  2) /6  =  669920. 


Figure  10.  A  Plot  of  Random  Data  and  Points  Forming  Three  Circular  Arcs 
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Figure  11.  Same  Plot  as  in  Figure  10  with  Circular  Arc  Points  Enhanced 
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The  extraction  of  the  plotted  points  of  the  three  L  — arcs  is  carried  out  by  the  subroutine 
CIRSOR  in  hie  DH.FOR.  In  more  detail,  the  generation  of  the  simulated  data  using  DOXY 
was  done  in  the  following  order:  First  N  (=160)  random  (x,y)  points  were  generated  and 
stored  in  array  AB;  then  alternating  starting  with  location  AB(5),  a  point  from  each  of  the 
L  — arcs  specihed  in  (31)- (33)  was  inserted  into  the  AB  overwriting  those  locations;  so  that 
L  —  al  is  stored  in  AB(5),  AB(8), . . . ,  AB(35);  L  —  a2  is  stored  in  AB(6),  AB(9), . . . ,  AB(66); 
La  —  3  is  stored  in  AB(7),  AB(IO), . . . ,  AB(52). 

The  output  form  CIRSOR  is  shown  in  Table  7. 

Table  7.  Listing  Output  Using  CIRSOR 

NR  =  112,  NC  =  48,  N  =  160.  See  (31)-(33)  for  L  — arc  points. 

TIME  29, 

L  -  al->  ICN(l)  =  11,  IC  =  5  8  11  14  17  20  23  26  29  32  35 

L  -  a2->  ICN(2)  =  21,  IC  =  6  9  12  15  18  21  24  27  30  33  36  39  42  45  48  51  54  57  60  63  66 

L  -  a3->  ICN(3)  =  16,  IC  =  7  10  13  16  19  22  25  28  31  34  37  40  43  46  49  52 


V.  EXTENSION  OF  THE  HOUGH  TRANSFORM  TO  ELLIPTICAL  ARCS 

We  consider  the  equation  for  the  ellipse  E  in  the  xy-plane  given  by 

(x-a)2  +  A(y-b)2  +  Bxy  +  D  =  0,  T  =  B^  -  4A  <  0.  (34) 

There  are  hve  unknowns  to  determine:  a,  b.  A,  B,  D.  We  use  hve  (x,  y)  points  onE  to  evaluate 
the  unknowns.  Or,  we  can  state  the  problem  as:  Given  hve  points  in  the  plane,  what  are  the 
values  of  the  unknowns  above  such  that  the  hve  points  are  onE?  Note  that  if  the  inequality 
in  (34)  is  not  satished,  then  other  geometrical  objects  will  be  determined  such  as  a  hyperbola, 
parabola,  or  straight  line.  However,  our  emphasis  will  be  with  ellipses. 

Let  the  hve  points  be  denoted  by  (xk,  yk),  k  =  1, ...,  5.  Then  four  linear  equations  can  be 
generated  for  a,  b.  A,  B  after  a^  and  b^  are  eliminated.  Namely 

2a(xk -xi)  +  A(yi  -  y^)  +  C  (yk -yi) +B  (xiyi -Xkyk)  =  Xk -xi,  k  =  2,  ...,5,  (35) 

where 

b  =  C/(2A).  (36) 

By  using  the  results  of  solving  the  linear  system  (35)  with  (36),  D  is  evaluated  from  (34) 
with,  say,  x  =  xi,  y  =  yi. 

Example  (calculations  carried  out  in  double  precision): 
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Specify  E'  by: 

50 

(x  -  5)2  +  — (y  -  5)2  +  2.5xy  -  94  |  =  0,  T  =  B2  -  4A  =  6.25  -  100/9  <  0.  (37) 

Thus, 

a  =  5,  b  =  5,  A  =  50/18,  B  =  2.5,  D  =  -94  |.  (38) 

The  objective  is  to  get  good  numerical  approximations  to  these  quantities  by  assigning 
hve  points  on  E.  Hence,  using  the  hve  points  onE,  (0,0),  (2,8.85078940),  (4,7.54511220), 
(6,6.03229152),  (8,4.17848880),  we  obtain  from  (35)  the  linear  system: 


/  2.00000000  -78.3364731  8.85078940  -17.7015788  \ 

/  2a  \ 

(  ^  ^ 

4.00000000  -56.9287181  7.54511220  -30.1804488 

A 

16 

6.00000000  -36.3885410  6.03229152  -36.1937491 

C 

36 

8.00000000  -17.4597686  4.17848880  -33.4279104  yi 

V  B  ) 

\qa) 

The  solution  of  (39)  with  (36),  as  obtained  by  subroutine  ELLP  in  ELLP5.FOR,  agrees  with 
the  values  of  those  quantities  as  given  in  (38)  to  within  eight  signihcant  digits.  Finally,  using 
those  results  in  (34),  with  xi  =  0,  yi  =  0,  gives  D  with  comparable  accuracy.  The  center 
point  coordinates  are  given  by  (xc,yc)  =  (2A[Bb  —  2a]/T,  2[a  —  2Ab]/T)  =  (—20/7,44/7). 
A  plot  oiE  is  shown  in  Figure  12  with  its  center  and  with  the  hve  specihed  points. 


Figure  12.  Plot  of  E  with  Five  Specihed  Points  in  xyp 
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VI.  FINDING  COORDINATES  OF  A  TARGET  FROM  TDOAs 


The  TDOAs  are  used  to  determine  the  coordinates  (Tx,Ty)  of  a  target  point  T,  by  using 
a  pulse  emitted  by  a  transmitter. 

Given  a  transmitter  Tl,  located  at  (Tlx,  Tly)  and  three  sensors,  Slk,  k=l,2,3  with 
locations  (Slkx,  Slky).  See  Figure  13. 


Figure  13.  Multistatic  Network  with  One  Transmitter  and  Three  Sensors 


A  pulse  PI  is  sent  from  Tl  to  T  and  received  from  T  at  Sll,  S12,  S13.  Let  the  time  taken 
for  each  of  these  paths  be  denoted  by  TIT,  TSll,  TS12,  and  TS13,  respectively.  Then  we 
can  dehne  times  of  arrival  TOAll  and  TOA12  by 

TOAll  =  TIT  +  TSll,  (40) 

TOA12  =  T1T  +  TS12.  (41) 

Subtracting  (40)  from  (41)  dehnes  TDOAl, 

TDOAl  =  TOA12  -  TOAll  =  TS12  -  TSll.  (42) 

Also, 

TOA13  =  TIT  +  TS13.  (43) 

Then,  subtracting  (40)  from  (43)  dehnes  TDOA2  accordingly 

TDOA2  =  TOA13  -  TOAll  =  TS13  -  TSll.  (44) 

We  are  now  in  position  to  achieve  the  objective  of  determining  the  target  point  coordinates 
(Tx,Ty).  Let  c  denote  the  speed  of  sound  in  meters/second,^  and  note  that 

c  *  TSlk  =  V(Tx  -  Slkx)2  +  (Ty  -  Slky)2,  k  =  l,2,3.  (45) 

=  331.3  *  +  X/273.15,  X  =  5/9  *  (F  —  32),  F  =  Degrees  Fahrenheit. 
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Then  from  (42)  and  (44)  we  have  two  nonlinear  equations  to  solve  for  (Tx,Ty),  namely 
c  *  TDOAl  =  v/(Tx  -  S12x)2  +  (Ty  -  S12y)2  -  ^(Tx  -  Sllx)2  +  (Ty  -  Slly)2,  (46) 

c  *  TD0A2  =  v/(Thr^Sm)^TTT^^^^Sl^  -  ^/{^x  -  Sllx)2  +  (Ty  -  Slly)2.  (47) 

Note  that  these  equations  are  independent  of  the  location  of  Tl. 

In  Fortran  hie  EXLAST2.F0R,  the  subroutine  HERD  [5]  is  used  to  solve  (46)  and  (47) 
for  Tx  and  Ty.  A  numerical  example,  assuming  acoustical  instruments,  with  distances  in 
meters  and  velocity  in  meters/seconds,  follows.  Input  is: 

TEMP  =  70°F,  c  =  343.8644, 

Sllx  =  3D3,  Slly  =  0,  S12x  =  5.71D3,  S12y  =  3.525D3,  S13x  =  1D3,  S13y  =  4D3,^ 
TDOAl  =  ID  -  1,  TD0A2  =  1.52D  -  1. 

With  an  initial  guess  of  Tx  =  Ty  =  1D3,  one  obtains,  as  output,  the  hnal  target  coordinates: 
Tx  =  3.2466935D3,  Ty  =  2.5890217D3. 

A  second  numerical  example  is  based  on  Jack  Carr’s  experimental  setup,  where  a  con¬ 
straint  is  imposed,  i.e., 

Ty  =  a  *  Tx  +  b,  (a,  b  given).  (48) 

In  this  case  only  two  sensors  are  required  resulting  in  one  TDOA. 

TEMP  =  70°F,  c  =  343.8644,  a  =  b  =  1, 

Sllx  =  3D3,  Slly  =  0,  S12x  =  5.71D3,  S12y  =  3.525D3, 

TDOAl  =  8D  -  2. 

Using  Fortran  subroutine  HERD  in  EXCONST. FOR,  with  an  initial  guess  of  Tx  =  1D3,  one 
obtains  the  hnal  target  coordinates: 

Tx  =  2.87598D3,  Ty  =  2.876.98D3. 

In  the  third  and  hnal  example,  RF  signals  are  transmitted  that  travel  at  the  speed  of 
light,  namely  c  =  2. 9979248D5  km/sec;  the  analysis  would  hold  for  acoustic  signals  as  well. 
Three  sensors  and  one  transmitter  make  up  the  network  to  locate  a  target  in  3-space  with  a 
constraint.  The  Earth  is  assumed  to  be  spherical  with  radius  Re  =  6378.137km. 

The  method  of  solution  to  determine  the  target  coordinates  (x,  y,  z)^  is  diherent  than  the 
one  used  in  the  previous  examples.  It  is  described  in  the  Ho-Chan  paper  [4]. 

Specihcally,  our  objective  is  to  develop  Fortran  95  software  that  is  based  on  one  of  the 
results  of  the  Ho-Chan  paper,  [4].  The  software  will  output  the  coordinates  of  a  target  from 
the  use  of  three  sensors,  whose  three-dimensional  position  coordinates  are  known,  with  the 
constraint  that  the  normal  distance  from  the  center  of  the  Earth  to  the  target,  r,  is  also 

^Numerical  values  are  specified  in  Fortran  notation.  For  example,  2. 3D  —  3  =  2.3  *  10“^  =  .0023.  All  computations  are 
carried  out  in  double  precision. 

^Hereafter  the  notation  of  [4]  is  used,  including  bold  upper  case  letters  to  denote  matrices  and  bold  lower  case  letters  to 
denote  vectors. 
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known.  As  noted  in  [4],  such  a  constraint  has  practical  application  as,  for  example,  a  target 
sitting  on  the  Earth’s  surface.  The  constraint  allows  the  location  of  the  target  with  only  two 
measurements  of  TDOAs. 

The  location  of  the  target  is  denoted  by  u  =  [x,  y,z]"'"  and  the  locations  of  the  three 
sensors  are  specihed  by  Sk  =  [xk,  yk,  Zk]"*",  k  =  1,2,3.  Thus,  the  distance  between  target 
and  sensor  k  is: 

rk  =  |sk  -  u|  =  (xk  -  x)2  +  (yk  -  y)^  +  (zk  -  z)^,  k=  1,2,3.  (49) 

To  clarify  the  meaning  of  TDOA  here,  consider  an  Earth  station  from  which  a  signal,  at 
arbitrary  time=0,  moving  at  the  speed  of  light,  c,  travels  to  the  target  and  the  reflected  signal 
then  travels  to  each  of  the  three  sensors  (It  is  assumed  that  not  any  three  of  {0,81,82,83} 
he  on  a  straight  line;  the  Earth’s  center  is  at  0.).  Let  TOAk  denote  the  time  of  arrival  of  the 
signal  at  sensor  k.  Then  let 

dk,i  =  TDOAk,i  =  TOAk  —  TOAi,  k  =  2,  3,  so  that 

rk,i  =  cdk,i  =  rk  -  ri,  k  =  2,  3.  (50) 

Note  that  the  constraint  can  also  be  written  as 


u^u  =  (Re  +  D)^  =  r^ 


(51) 


where  Re  =  spherical  Earth's  radius,  D  =  distance  above  the  Earth  along  u. 

From  this  point,  we  proceed  mathematically,  following  [4],  toward  the  objective  of  Ending 
u,  given:  8k,  TDOA2,i,  TDOAa^i,  and  r.  Rewriting  (50)  as  rk,i  +  ri  =  rk  and  squaring  both 
sides  of  this  relation  gives 


‘■k.i 


2rk,iri  +  ri  =rk  =  (sk-u)'^(8k-u)  =r^  +  8kSk-28kU,  k  =  2,3,  (52) 


where  ri  from  (49)  gives 


ri  =  r 


87  81  —  287  u. 


Replacing  on  the  left  side  of  (52)  by  (53)  gives 


k,l 


2rk,iri  =  sjsk  -  sj'si  -  2(sk  -  Si)’'u,  k  =  2,3, 


(53) 


(54) 


Equations  (53)  and  (54)  represent  a  set  of  linear  equations  for  u  in  terms  of  the  variable  ri. 
In  matrix  notation  we  have 

Gu  =  h,  (55) 

where 


G  =  -2 


T  T 

8}  -  8} 

T  T 

S3  - 


(56) 


20 


NSWCDD/TR-08/136 


1  -r^  -  s^Si  Oil 

h  =  H  n  =  -  sjs2  +  s^Si  2r2,i  0  ri  .  (57) 

jl_  _r|  i  -  sjs3  +  s^si  2r3,i  0_  jf_ 

As  noted  in  [4],  the  existence  of  requires  that  any  three  of  (0,  Si,  S2,  S3)  cannot  he  on  a 
straight  line.  In  addition,  the  three  sensors  should  have  sufficient  spatial  separation  to  avoid 
ill-conditioning  of  G.  Hence,  from  (55), 


where 


Substituting  into  (51),  we  have 


u  =  G  =  Q  [l,ri,ri 


Q  =  G 


T  2 

u  u  —  r 


=  [l,ri,rnQ^Q[l,ri,r?r-r2  =  0. 


which  represents  a  quartic  polynomial  in  ri.  Let  P  denote  a  symmetric  matrix  such  that 

Pii  P12  Pl3 

P  =  P12  P22  P23  =  Q"'"Q  (61) 

Pl3  P23  P33 

Expanding  (60)  and  using  (57),  (59),  and  (61),  we  obtain  the  coefficients  of  the  quartic  F. 
Thus, 

F(ri)  =  ( pii  -  r^  )  +  2  P12  ri  +  ( 2  P13  +  P22  )rl  +  2  P23  rf  +  P33  r^  =  0  (62) 

We  are  interested  in  the  positive  roots  of  F.  Using  one  such  root  in  (58)  gives  a  set  of 
possible  target  coordinates.  In  case  there  is  more  than  one  positive  root,  as  pointed  out  in 
[4],  extraneous  roots  will  not  satisfy  (50).  However,  it  still  can  happen  that  more  than  one 
positive  root  satishes  (50),  in  which  case  the  user  must  decide  by  other  means  which  is  the 
meaningful  one  for  his  application.  The  numerical  example  that  follows  reflects  such  a  case. 

The  hie  CHAN. FOR  uses  four  subroutines  from  the  NSWC  Library  of  Mathematical  Sub¬ 
routines,  [5]: 

(1)  CROUT  (Input:  G,  output:  G^^  stored  in  G.) 

(2)  MTMS  (Multiplies  G-^  H,  stores  result  in  Q.) 

(3)  TMPROD  (Input:  Q,  output:  Q^Q  stored  in  QT.) 

(4)  DRPOLY  (Input:  Coefficients  of  polynomial  F  (elements  of  QT,  see  (62))  stored  in  Tl, 
output:  real  parts  of  roots  of  F  in  array  ZR  and  imaginary  parts  in  array  ZI.) 

We  use  c  (speed  of  light)=  2.9979248D5  km/sec.  Re  (radius  of  the  Earth  )=  6378.137  km. 
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Name  of  Fortran  95  file:  Chan. for  (Contains  main  program  plus  the  four 
mentioned  above,  including  their  supporting  routines.) 

Input:  D  (distance  above  the  Earth  of  the  target)=  1  km 


Si 

=  [0.0 

0.0 

6378. 137]^ 

S2 

=  [3189.069 

3189.069 

4510.024]^ 

S3 

=  [2761.814 

4783.603 

3189.069]'^ 

TDOAi,2 

=  2D -4 

TDOAi,3 

=  5D-4 

r 

=  6379.137 

Intermediate  Results: 


G 


H 


7.81071D  -  13 

-6.37814D3 

-5.52363D3 


O.OOOOODO 

-6.37814D3 

-9.56721D3 


-1.27563D4 

3.73623D3 

6.37814D3 


-8.13740D7  O.OOOOODO  l.OOOOODO 
3.59502D3  1.19917D2  O.OOOOODO 

2.24689D4  2.99792D2  O.OOOOODO 


G  1 


1.50016D  -  5 
-6.09230D  -  5 
-7.83928D  -  5 


-3.70959D  -  4 
2.14173D  -4 
2.27139D  -  20 


2.47306D  -  4 
-2.47306D  -  4 
-1.51426D  -  20 


Q  = 


Q^Q  = 


-1.21652D3  2.96562D  -  2  1.50016D  -  5 

4.95277D3  -4.84574D  -  2  -6.09230D  -  5 

6.37914D3  -1.81586D  -  18  -7.83928D  -  5 

6.67032D7  -2.76076D2  -8.20066D  -  1 

-2.76076D2  3.22761D  -  3  3.39706D  -  6 

-8.20066D  -  1  3.39706D  -  6  1.00821D  -  8 

Det  =  — 3.28991D11  (Determinant  of  G) 


where 


F(ri)  =  5^Tl(k)r^"  =  0, 

k=i 

Tl(l)  =  1.00821D-8 
Tl(2)  =  6.79413D-6 
Tl(3)  =  -1.63690D0 
Tl(4)  =  -5.52151D2 
Tl(5)  =  2.60098D7. 


subroutines 


(63) 


(64) 


(65) 


(66) 
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With  F  as  given  in  (65),  DRPLOY  finds  the  roots: 


ZR(1)  = 

4.05816D3, 

ZI(1) 

=  0.0 

ZR(2)  = 

-4.39528D3, 

ZI(2) 

=  0.0 

ZR(3)  = 

1.18592D4, 

ZI(3) 

=  0.0 

ZR(4)  = 

-1.21960D4, 

ZI(4) 

=  0.0 

Output: 

We  are  only  interested  in  the  positive  roots  ZR(1)  and  ZR(3).  Noting  from  (58)  that 

u  =  Q[l,ri,r2]^ 

there  are  two  solutions,  namely 

u  =  [-8.49113D2,  3.75280D3,  5.08811D3]^ 

u  =  [1.24501D3,  -4.19015D3,  -4.64607D3]'^. 

Both  solutions,  using  (49),  satisfy  (50)  with 

r2,i  =  59.95848,  rg^i  =  149.8962. 

Thus,  neither  positive  root  of  (65)  is  extraneous  and,  therefore,  both  solutions  hold. 
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APPENDIX  A 

LIST  OF  FORTRAN  95  SUBPROGRAMS 

1.  MCARR.FOR  or  MCARRl. FOR— File  containing  programs  to  generate  L  — lines  from 
an  M-set,  based  on  the  Hough  Transformation. 

(a)  GUVM— Subroutine  that  generates  an  M-set  of  points  to  test  HSORTl. 

(b) HSORTl— Subroutine  that  Ends  L  — lines  in  an  M-set  of  points. 

(c) MATPLO— Generates  a  Matlab  M  hie  for  plotting  results  from  HSORTl. 

(d) HSORTI— Sorting  routine  from  Mathlib,  [5]  of  the  main  text. 

2.  DH. FOR— File  containing  subprograms  to  generate  L  — arcs  from  an  M-set. 

(a)  DGXY— Subroutine  generating  an  M-set  of  points  for  testing  subroutine  GIRSOR. 

(b) GIRl— Subroutine  that  finds  parameters  defining  a  circle  from  three  M-points. 

(c)  GIRSOR— Subroutine  that  finds  L  — arcs  from  an  M-set  of  points. 

(d) MATPLX— Generates  an  Matlab  M  hie  for  plotting  results  from  GIRSOR. 

(e) LE— Supporting  routine  for  GIRSOR.  See  Appendix  B. 

(f) LGALl— Supporting  routine  for  LE. 

(g) RNORM— Normal  random  number  generating  routine  used  in  DGXY,  taken  from 
[5]  of  the  main  text. 

(h) RGIRl— Supporting  routine  for  RNORM. 

3.  PTGIRG. FOR— File  containing  subprograms  to  generate  L— lines,  each  with  an  attached 
circular  arc  at  one  endpoint. 

(a)  GXYG— Generates  an  M-set  of  points  to  test  HSORTl  together  with  HS0RT2. 

(b) HSORTl— See  above. 

(c) HS0RT2— From  an  M-set,  finds  a  circular  arc  attached  to  an  endpoint  of  an  L— line. 

(d)  GUV— Subroutine  that  generates  the  parameters  of  an  L  — line  of  an  M-set. 

(e) GUVR— Subroutine  that  generates  the  parameters  of  a  circular  arc  attached  to  an 
L  — line. 

(f) GIR— Subroutine,  essentially  the  same  as  GIRl  above. 

(g) MATPLO— Generates  a  Matlab  M  hie  for  plotting  results  from  GXYG. 

(h) RNORM— See  above. 

(i) RGIRl— See  above. 

4.  ELLP5. FOR— File  containing  subprograms  to  find  parameters  defining  an  ellipse,  given 
five  points  in  the  xy-plane. 

(a)ELLPT— Given  ellipse  parameters  and  ellipse  coordinates  xl,yl,x2,  ELLPT  finds 
two  y-coordinates  on  the  ellipse  with  x2  as  the  x-coordinate. 
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(b) ELLP— Given  five  points  in  the  plane,  ELLP  determines  the  ellipse,  E,  with  those 
points  on  E,  by  outputting  the  parameters  dehning  E. 

(c) CROUT— Subroutine  that  solves  a  linear  system,  with  real  coefficients,  by  the  Grout 
procedure.  Gontained  in  [5]  of  the  main  text. 

(d)  ROTA— Gives  the  angle  of  rotation  of  the  xy-axes  so  that  the  xy  term  in  (34)  of 
the  main  text  is  eliminated.  The  output  also  includes  the  parameters  defining  E  in  the 
rotated  coordinates. 

(e) RNORM— See  above. 

(f) RGIRl— See  above. 

5.  EXLAST2. FOR— File  containing  subprograms  for  determining  target  coordinates  from 
a  specified  set  of  sensors  and  a  transmitter. 

(a)HBRD— Subroutine  that  solves  a  set  of  nonlinear  functions.  Gontained  in  [5]  of  the 
main  text. 

6.  GHAN. FOR— File  containing  subprograms  for  using  the  Ghan  method  to  determine 
three-dimensional  target  coordinates  using  three  sensors  and  a  constraint. 

(a) GROUT— See  above. 

(b) MTMS— Subroutine,  from  [5]  of  the  main  text,  for  multiplying  two  real  matrices. 

(c) TMPROD— Subroutine,  from  [5]  of  the  main  text,  to  produce  A"*^  B  given  real  ma¬ 
trices  A  and  B. 

(d) DRPOLY— Subroutine,  from  [5]  of  the  main  text,  to  produce  the  roots  of  a  polyno¬ 
mial  with  real  coefficients. 
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APPENDIX  B 

A  RESULT  USED  IN  SUBROUTINE  LCALl  CALLED  BY  SUBROUTINE  LE 

This  analysis  is  used  in  subroutine  LE  of  file  HCIR.FOR  and  the  double  precision  version 
DHCIR.FOR.  It  is  also  used  in  the  latest  circular  arc  routine  DCIR3.FOR.  LCALl  replaces 
the  subroutine  IJKL,  which  is  much  slower. 

We  have  the  loops: 

L  =  0 

DO  5  1=1,  N-2  !N  IS  TOTAL  NO.  OF  (X,Y)  POINTS.  SEE  ROUTINE  CIRSOR 
DO  5  J=I+1,N-1 
DO  5  K=J+1,N 
L  =  L  +  1 
5  CONTINUE 

PROBLEM:  Given  N,  I,  J,  K.  Find  L  without  running  through  3  LOOPS. 

For  example:  N=10, 1=2,  J=5,  K=6;  FIND  L.  It  is  L=50. 

We  use: 

n  n 

m  =  n(n  +  l)/2,  m^  =  n(n  +  l)(2n  +  l)/6.  (B-1) 

m=l  m=l 

We  first  hud  the  contribution  to  L  from  I;  call  it  Lj.  Then 


1  ^ 

Lj  = -X(N-m)(N-m- 1),  1  =  1-1. 

(B-2) 

2  Lj 

m — 1 

=  i  n2  -  i  N  (i  + 1)  +  i  i  (i  + 1)(2  i  + 1)  -  i  N  +  i  i  (i  + 1), 

0  2 

(B-3) 

2  Lj 

=  iN^  -  IN  (i  +  2)  +  i  i  (i  +  i)(2i  + 1)  +  ^  i  (i  + 1), 

D  2 

(B-4) 

2  Lj 

=  iN(N-i-i)  +  lii(i  +  i), 

(B-6) 

^  i  [3N(N-I-1)  +  I(I  +  1)] 

Lt  — 

^  6 

(B-6) 

For  a  fixed  I,  the  contribution  to  L  from  J  and  K,  denoted  by  Ljk,  is  given  by 

j-i 

Ljk  =  XI  K  -  J’ 

n=I+l 

Ljk  =  (J-I-l)N  -  J(J-l)/2  +  I(I+l)/2  +  K- J, 

Ljk  =  [(J-I-1)(2N-I- J)]/2  +  K- J. 
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Hence 

L  =  Li  +  Ljk-  (B-10) 

Completing  the  numerical  example  above,  we  get  for  N  =  10,  I  =  2,  J  =  5,  K  =  6; 

Li  =  36,  Ljk  =  14,  L  =  50. 
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